//	--------------------------------------------------------------------
//	Copyright(C) 2006,2007 Zhao Yukun. All rights reserved.
//	This file is a part of the Xeres system.
//	Xeres 2007
//	Contact : xeres.engine@gmail.com
//	$(reservedInfo)
//	$(reservedDeclaration)
//	--------------------------------------------------------------------
///	@file	<random.h>
///	@path	~/src/lib/math/
///	@date	2008/05/15
///	@desc	Random number generator.

#pragma once

#include "xeres.h"

namespace xeres {

	/* unit random number generator [0,1] ( depend on rand() )*/
	inline const float urand( void )
	{
		static const float r = (float)( 1. / (double)RAND_MAX );
		return (float)rand() * r;
	}

	/* ---------------------------------------------------
	
		Mersenne Twister Pseudo Random Number Generator
	
	--------------------------------------------------- */
	class MTRandom
	{
	private:
		static const int mt19937_n						= 624;
		static const int mt19937_m						= 397;
		static const unsigned int mt19937_matrix_a		= 0x9908b0dfUL;
		static const unsigned int mt19937_upper_mask	= 0x80000000UL;
		static const unsigned int mt19937_lower_mask	= 0x7fffffffUL;
	private:
		unsigned int mt[mt19937_n]; /* the array for the state vector  */
		int mti; /* mti==mt19937_n+1 means mt[mt19937_n] is not initialized */
	public:
		MTRandom( unsigned int seed = 0 ) {
			reseed( seed );
		}
	public:
		/* initializes mt[mt19937_n] with a seed */
		void reseed( unsigned int seed ) {
			mt[0]= seed & 0xffffffffUL;
			for( mti = 1 ; mti < mt19937_n ; mti++ )
			{
				mt[mti] = ( 1812433253UL * ( mt[mti - 1] ^ ( mt[mti - 1] >> 30 ) ) + mti );
				/* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
				/* In the previous versions, MSBs of the seed affect   */
				/* only MSBs of the array mt[].                        */
				/* 2002/01/09 modified by Makoto Matsumoto             */
				mt[mti] &= 0xffffffffUL;
				/* for >32 bit machines */
			}
		}
		/* initialize by an array with array-length */
		/* init_key is the array for initializing keys */
		/* key_length is its length */
		void reseed( unsigned int init_key[] , size_t key_length ) {
			size_t i, j, k;
			reseed( 19650218UL );
			i = 1; j = 0;
			k = ( mt19937_n > key_length ? mt19937_n : key_length );
			for( ; k ; k--)
			{
				mt[i] = ( mt[i] ^ ( ( mt[i-1] ^ ( mt[i - 1] >> 30 ) ) * 1664525UL ) ) + init_key[j] + j; /* non linear */
				mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
				i++; j++;
				if( i >= mt19937_n ) { mt[0] = mt[mt19937_n - 1]; i = 1; }
				if( j >= key_length ) j = 0;
			}
			for( k = mt19937_n - 1 ; k ; k-- )
			{
				mt[i] = ( mt[i] ^ ( ( mt[i - 1] ^ ( mt[i - 1] >> 30 ) ) * 1566083941UL ) ) - i; /* non linear */
				mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
				i++;
				if( i >= mt19937_n ) { mt[0] = mt[mt19937_n - 1]; i = 1; }
			}
			mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */ 
		}
	public:
		/* generates a random number on [0,0xffffffff]-interval */
		unsigned int randu32( void ) {
			unsigned int y;
			static unsigned int mag01[2] = { 0x0UL , mt19937_matrix_a };
			/* mag01[x] = x * MATRIX_A  for x=0,1 */
			if( mti >= mt19937_n )
			{ /* generate mt19937_n words at one time */
				int kk;

				if( mti == mt19937_n + 1 )   /* if init_genrand() has not been called, */
					reseed( 5489UL ); /* a default initial seed is used */
				for( kk = 0 ; kk < mt19937_n - mt19937_m ; kk++ )
				{
					y = ( mt[kk] & mt19937_upper_mask ) | ( mt[kk + 1] & mt19937_lower_mask );
					mt[kk] = mt[kk + mt19937_m] ^ (y >> 1) ^ mag01[y & 0x1UL];
				}
				for( ; kk < mt19937_n - 1 ; kk++ )
				{
					y = ( mt[kk] & mt19937_upper_mask ) | ( mt[kk + 1] & mt19937_lower_mask );
					mt[kk] = mt[kk + ( mt19937_m - mt19937_n )] ^ (y >> 1) ^ mag01[y & 0x1UL];
				}
				y = ( mt[mt19937_n - 1] & mt19937_upper_mask ) | ( mt[0] & mt19937_lower_mask );
				mt[mt19937_n - 1] = mt[mt19937_m - 1] ^ (y >> 1) ^ mag01[y & 0x1UL];
				mti = 0;
			}

			y = mt[mti++];
			/* Tempering */
			y ^= ( y >> 11 );
			y ^= ( y << 7 ) & 0x9d2c5680UL;
			y ^= ( y << 15 ) & 0xefc60000UL;
			y ^= ( y >> 18 );
			return y;
		}
		/* generates a random number on [0,0x7fffffff]-interval */
		int randi31( void ) {
			return ( randu32() >> 1 );
		}
		/* generates a random number on [0,1]-real-interval */
		double randd1( void ) {
			static const double f = 1.0/4294967295.0;
			return ((double)randu32()) * f;
			/* divided by 2^32-1 */ 
		}
		/* generates a random number on [0,1)-real-interval */
		double randd2( void ) {
			static const double f = 1.0/4294967296.0;
			return ((double)randu32()) * f;
		}
		/* generates a random number on (0,1)-real-interval */
		double randd3( void ) {
			static const double f = 1.0/4294967296.0;
			return ((double)randu32() + 0.5 ) * f;
		}
		/* generates a random number on [0,1) with 53-bit resolution*/
		double randd53( void ) {
			static const double f = 1.0/9007199254740992.0;
			unsigned int a = randu32() >> 5 , b = randu32() >> 6;
			return ( (double) a * 67108864.0 + b ) * f; 
		}
		/* generates a random number on [0,1]-real-interval */
		float randf1( void ) {
			static const float f = (float)(1.0/4294967295.0);
			return (float)((float)randu32()) * f;
			/* divided by 2^32-1 */ 
		}
		/* generates a random number on [0,1)-real-interval */
		float randf2( void ) {
			static const float f = (float)(1.0/4294967296.0);
			return (float)((float)randu32()) * f;
		}
		/* generates a random number on (0,1)-real-interval */
		float randf3( void ) {
			static const float f = (float)(1.0/4294967296.0);
			return (float)((float)randu32() + 0.5 ) * f;
		}
	};

} // namespace xeres
